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We discuss spin-| one-dimensional (ID) and quasi-lD magnets with competing ferromagnetic 
nearest-neighbor Ji and antiferromagnetic next-nearest-neighbor J exchange interactions in a strong 
magnetic field H. It is well known that due to attraction between magnons quantum phase transi- 
tions (QPTs) take place at H = H s from the fully polarized phase to nematic ones if J > |Ji|/4. 
Such a transition at J > 0.368| J\\ is characterized by a softening of the two-magnon bound-state 
spectrum. Using a bond operator formalism we propose a bosonic representation of the spin Hamil- 
tonian containing, aside from bosons describing one-magnon spin-1 excitations, a boson describing 
spin-2 excitations whose spectrum coincides at H > H 3 with the two-magnon bound-state spectrum 
obtained before. The presence of the bosonic mode in the theory that softens at H — H s makes 
the QPT consideration substantially standard. In the ID case at H < H s , we find an expression 
for the magnetization which describes well existing numerical data. Expressions for spin correlators 
are obtained which coincide with those derived before either in the limiting case of J > |Ji| or 
using a phenomenological theory. In quasi-lD magnets, we find that the boundary in the H—T 
plane between the fully polarized and the nematic phases is given by H a (0) — H a (T) oc J 13 / 2 . Sim- 
ple expressions are obtained in the nematic phase for static spin correlators, spectra of magnons 
and the soft mode, magnetization and the nematic order parameter. All static two-spin correlation 
functions are short ranged with the correlation length proportional to l/ln(l + | J\\/ J). Dynamical 
spin susceptibilities are discussed and it is shown that the soft mode can be observed experimentally 
in the longitudinal channel. 

PACS numbers: 75.10.Jm, 75.10.Kt, 75.10.Pq 
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FIG. 1. (Color online.) (a) Quasi-ID frustrated magnet described by Hamiltonian |T]). It is implied that inter-chain interactions 
J[ and J2 are small compared to those inside chains. The approach suggested in the present paper is based on the unit cell 
doubling along chains direction that is shown in panel (b). 

I. INTRODUCTION 

Spin nematic states with multiple-spin ordering and without the conventional long-range magnetic order have 
attracted much attention in recent years. Two-spin nematic order can be generally described by the tensor^ = 

(S"Si) — 5 a p(SjSi) /3. The antisymmetric part of is related to the vector chirality (Sj x S/). The formation of 
the vector chiral spin liquid was anticipated a long time ago in 2D frustrated spin systemsP Such states have been 
obtained recently in a ring-exchange spin-| model at T = 0® and in classical frustrated spin systems at T / OP 
The symmetric part of describes a quadrupolar order. In this instance one distinguishes the cases of j = I and 
j =/= I. It has been known for a long time that the one-site (j — I) nematic state can be stabilized by the sufficiently 
strong biquadratic exchange (SiS2) 2 .'^The interest in this mechanism of multiple-spin ordering stabilization has been 
revived recently in connection with experiments on cold atom gases^and on the disordered spin system NiGa2S4pl 

While the one-site quadrupolar ordering can exist only for S > 1, the different-site one (j ^ I) can be found even 
in spin- \ systems. In particular, the different-site nematic phases have been discussed recently in (quasi-) ID 
2DJ22HH] anc | 3E(2Sl systems with competing ferro- and antiferromagnetic exchange couplings between neighboring and 
next-neighboring spins, respectively. The (quasi-) ID magnet of this kind is described by the Hamiltonian 

u = (-SA+i + JSjS j+2 ) -hJ2s] + h', (1) 

where we set the ferromagnetic exchange coupling constant between neighboring spins to be equal to —1 and H' 
describes an inter-chain interaction that is also taken into account in the present paper (see Fig. [lja)). As the field 
direction can be arbitrary, we direct the field perpendicular to chains for simplicity. The interest in model ([!]) is 
stimulated also by recent experiments on the corresponding quasi-lD materials LiCuVO^J^HaU Rb 2 Cu2Mo30i2!^ 
Li 2 ZrCu0 4 , [SI1 CuCl2, ISS1 PbCuS04(OH)2, l33 LiCuSb04, ISS1 and some others. 

It has been found recently that the physics of spin-^ model with %' = is even richer: field-driven transitions 
to phases with quasi-long-range multiple-spin ordering have been obtained below the saturation field H s and are 
described by operators 5 - S- +1 . . . Sf +p _ 1 with p = 2 (quadrupolar phase) at J > 0.368, p = 3 (hexapolar phase) at 
0.284 < J < 0.368 and p = 4 (octupolar phase) at 0.259 < J < 0.284PM 1 Finite inter-chain interaction W stabilizes 
the long-range nematic order at T = but quite a small nonfrustrated T~i! turns the point H = H s in to an ordin ary 
quantum critical point separating the fully polarized phase and that with a long-range magnetic order jST 18 ! 21 ! 22 ! 

It is well known that the origin of the nematic phases is the attraction between magnons caused by frustration.^ 
As a result of this attraction, the bottom of the one-magnon band lies above the lowest multi-magnon bound state 
at H = H s (see, e.g., Fig. [2]ja)). Then, transitions to nematic phases are characterized by a softening of the multi- 
magnon bound-state spectrum rather than the one-magnon spectrum. As a consequence, new approaches are required 
to describe such transitions. 
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FIG. 2. (Color online.) (a) Spectra of one-magnon excitations and two-magnon bound states are shown at H — H s as they 
were found before at J > 0.375 for an isolated chain described by Hamiltonian |l| (the particular curves drawn are for J = 1). 
The distance between two neighboring spins in the chain is equal to d. Dashed lines are for parts of the bound-state spectrum 
which lie in the two-magnon continuum. The first Brillouin zone (BZ) [— n/2d, n/2d] is marked after the unit cell doubling that 
is drawn in Fig. [ljb). The result of the reduction to this BZ is presented in panel (b). Within the approach suggested in the 
present paper, spectra of spin-1 excitations described by a- and ft- particles (see Eqs. |6|-( 12 |) are two parts of the one-magnon 



spectrum shown in panel (a), whereas the spectrum of a-particles (spin-2 excitations) coincides with the low-energy part of the 
two-magnon bound-state spectrum. 



Properties of model ([I]) are well-understood in the fully polarized state at H > H s . Wave functions of the multi- 
magnon bound states can be represented as linear combinations of functions S~S~ . . . |0), where |0) is the vacuum 
state at which all spins have the m aximu m projection on the field direction, and the spectrum can be found numerically 
from the corresponding equations! 16 * 39 * An approach is suggested in RefP^ that allows one to map the system at 
H > H s to a tight-binding impurity problem. In particular, it allows authors to obtain analytical expressions for 
the two-magnon bound-state spectrum and to show that it is quadratic at J > 0.375 near its minimum located at 
k = ir/d, where d is the distance between neighboring spins (see Fig. Kla)). 

Transitions to nematic spin liquid phases in the purely ID spin-~| model are discussed using the bosonization 
technique in the limit J ^ 1 and using a phenomenological approach at arbitrary J > 1 /4.HHIir i3 l 16 l According to the 
latter method the transition is equivalent to that in ID hard-core Bose-gas with the following correspondence between 

the bosonic operators b, and spins: b\ = (-l) l S^S~ +1 . . . S'^ f) _ 1 = (-l) l Aff p) and b\bi = (1/2 - S?)/p. Results of 
the phenomenological approach agree with those of the bosonization method in the region of its validity J>1 and 

show, in particular, an algebraic decay of static spin correlators (S'fS'J) and (M^Mj } in nematic phases and an 
exponential decay of (S^~ S~). Although many predictions of the phenomenological theory are confirmed by numerical 
calculations, the corresponding microscopic analytical calculations based on the spin Hamiltonian are also desirable 
at J ~ 1. 

It is well known that the behavior of quasi-lD systems differs significantly at low T from that of purely ID systems. 
Then, a special approach for a quasi-lD model Q at T = has been suggested recent lypH The wave function of the 
ground state in the quadrupolar phase is proposed in a form that resembles the BCS pairing wave function of electrons 
in superconductors. Using this approach authors have calculated static spin correlators and found that in contrast to 
the purely ID case {SfS?) decays exponentially and the system has a long-range nematic " antiferromagnetic" order. 
The magnetization of LiCuV04 at H < H s measured in Ref is also described successfully in RefPS. However many 
dynamical properties as well as the temperature effect have not been considered yet in the quasi- ID case that leaves 
room for further theoretical discussion in this field. 

We suggest an approach in the present paper that allows us to perform a quantitative microscopic consideration 
of the quadrupolar phase at T > and H s» H s both in purely ID and quasi- ID spin-| model ([!]). This approach 
is based on the unit cell doubling along the chain direction that is shown in Fig. |ljb) and on a representation of 
two spin operators in each unit cell via three bosons. This representation resembles those proposed for dimer spin-| 
systemsP^H^ Two of these bosons describe one-magnon (spin-1) modes while the third one describes spin-2 excitations 
(see Fig. [2]) which are referred to as a-particles below. We demonstrate that the spectrum of a-partic les coincides at 
H > H s with the spectrum of two-magnon bound states calculated before by other methods! 8 * 17 * 18 ! Then, it is the 
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main advantage of our approach that it contains a boson whose spectrum becomes "soft" as a result of the transition 
to the nematic phase. This circumstance makes relatively simple and quite standard the quantitative discussion of 
the transition. 

It should be noted at once that the procedure of unit-cell doubling is arbitrary (there are two ways for neighboring 
spins to group into couples) and it breaks the initial translational symmetry. However, it does not play a role in our 
consideration because all the physical results are obtained in the present paper either exactly or using perturbation 
theories with "good" small parameters. As a consequence, the translational symmetry is restored in our results at 
H > H s , it turns out to be broken at H < H s in accordance with conclusions of previous considerations and all the 
physical results obtained at H < H s do not depend on the way of spins grouping into couples. 

By using our approach, we confirm the hypothesis proposed in RefJ^that the transition in an isolated chain to the 
quadrupolar phase is equivalent in many respects t o that in ID systems of hard-core bosons. We rederive many of 
the results for the isolated chain obtained in RefsPS ^ 13 * 16 ! As an extension of the previous discussion, we derive an 
expression for the magnetization that describes well available numerical data at H ss H s . 

In quasi-lD systems, we find that the boundary in the H—T plane between the fully polarized and the quadrupolar 
phases is given by H s (0) - H S (T) cx T 3 / 2 . Simple expressions are obtained for static spin correlators, spectra of 
magnons and the soft mode, magnetization and the nematic order para meter. We obtain an " antiferromagnetic" 
nematic long-range order along the chains in accordance with Refs.^H. All the static two-spin correlators decay 
exponentially with the correlation length proportional to l/ln(l + 1/J). This exponential decay results in broad 
peaks in the transverse structure factor with the period along chains equal to n/d rather than 2tt jd. Dynamical spin 
susceptibilities x Q ^(w, q) are discussed, where a,j3 = x, y, z. It is shown that Xzz( w , q) has sharp peaks at w equal to 
energies of a-particles. Thus, the soft mode can be observed experimentally in the longitudinal channel. There are 
sharp peaks at uj corresponding to energies of magnons in transverse components of the dynamical spin susceptibility 
(in accordance with predictions of RefPS). An application is discussed of the proposed theory to LiCuVO.4. Our 
results are in reasonable agreement with available experimental data for magnetization in this compound. 

The rest of the present paper is organized as follows. We describe in detail our approach in Sec. [TTJ. Properties of an 
isolated chain are discussed in Sec. |III| Quasi-lD systems at H > H s and H < H s are considered in Sees. |IV| and[V] 
respectively. Sec. |VI| contains a summary of our results and a conclusion. Some details of calculations and the model 
describing LiCuV04 are discussed in appendices. 



II. APPROACH 



A. Spin representation 

We start with an isolated chain at H > H s . To describe its properties in the fully polarized state and the quantum 
phase transition we suggest to double the unit cell as it is shown in Fig.[TJb). Then, there are two spins, S\j and S23, 
in j-th unit cell. To take into account all spin degrees of freedom in each unit cell we introduce three Bose-operators 
Oj, b^ and c| which create three spin states from the vacuum |0) as follows: 

10} = |tt) , 
oj|0) = III) , 

&}|0) = III) , (2) 
c}|0) = lit), 

where all spins have the maximum projection on the field direction at the state |0). One leads to the following spin 
representation via these Bose-operators: 

S h = c ] a i + b h 

S 2j= a j c j + b p ( 3 ) 
$2j =2 - a j a i ~ h \ h r 

It is easy to verify that Eqs. ^ reproduce spin commutation relations on the physical subspace (which consists of 
states with no more than one particle a, b or c in each unit cell) of the Hilbert space and Sy = S| ■ =3/4. In order 
to eliminate contributions to physical quantities from unphysical states one can introduce into Eqs. O) the projector 



s h = b j a j+ c J> 
s ij = a ] b J + c p 

Sij = - — djdj — Cj c j, 
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operator 1 — oJjCtj — bjbj — c -Cj or add to the Hamiltonian a term describing infinite repulsion between particles in 
each unit cell 

U J2{ a l a ) a J a i + 4 4 b J h i + 4 4 c i C J + 4 4 a i b J + 4 4 a J c i + 4 4 b i c i ) ' U 00 ■ ( 4 ) 

Both methods should lead to the same results at small T (see, e.g., RefP) and we choose the last one in the present 
paper. 

Representation Q-Q is an analog of the bond-operator representation suggested in RefP^and applied to systems 
with singlet (" dimerized" ) ground states. In principle, Eqs. Q-Q can be derived from RefP^ implying that \\\) is 
the ground state as it is done, e.g., in RefP^lfor the singlet ground state. 



B. Hamiltonian transformation 

Substituting Eqs. ([3| into Hamiltonian with W = and taking into account constraint Q one obtains 

H = Sa + H2+H 3 +H4, (5) 

where is a constant, 

n 2 =Y, ( e a o4«k + £ k (jibk + c£ck) + Bkc^k + B k 6 k c k ) , (6) 

k 

7^7 £ f"I ( elkl + e< *") 44^-3 - \ (e- lkl + e-' lk2 ) 4^-3 + \{Ji + .h)b\c\a^ 



n 3 = 



N 

k!+k2+k3=0 



\ (e tk2 + e^) a t 1 c_ 2 c_ 3 ~ \ (e- lk2 + e~^) a\b_ 2 b_ 3 + \(J 2 + J 3 )aJ&_ 2 c_ 3 ) , 
\ (( U + J i+3 ~ e 4(fel+fc3) ) a t 1 4a_ 3 a_ 4 + (V + i J 1+3 \ (&I&|6_ 3 &_ 4 + C t 1 c^_ 3 c_ 4 ) 

k 1 +k 2 +k 3 +k 4 =0 > V / 

+ (u+ J1+3 + \ji+4 - e-^+^A 44a_ 3 c_ 4 + (u + Ji +3 + i Ji+4 - e^+^A a t 1 6 2 a_ 3 6_ 4 
+ (U- e - l ^ +k ^) 6 t 1 c|6_ 3 c_ 4 - ^- i(fe2+fe3) a t lC ^- 3 6-4 - \er^ +k ^ a\b\a^c^ , 



(7) 



(8) 



k is the one-dimensional momentum here, Jk = 2 J cos fc, we set 2d = 1, e a0 = 2_ff + 1 — J , = H + 1 — ( J — J k )/2, 
Bk = —e %k l 2 cos |, TV is the number of unit cells (that is half the number of spins in the lattice), and we omit some 
indexes k in Eqs. Q and M. 
After the unitary transformation 

c k = ~e ifc / 4 (a k + /3 k ), fo k - ^e- lfe / 4 (/3 k - a k ) , (9) 

the bilinear part of the Hamiltonian ^ acquires the form 

Ui = ( e ao«k a k + ea(k)a k a k + e^(k)/9^ k ) , (10) 

k 

where 

fc 

e Q (k) = H + 1 - J + Jcosfc + cos -, (11) 

k 

ep{k) = H + 1 - J + Jcosfc - cos -. (12) 



Eqs. ( 11 ) and ( 12 1 represent two branches of the one-magnon spectrum which result from the unit-cell doubling as it 
is shown in Fig. |2j The lower branch, ep(\s), has minima at k = (±Q, 0,0), where 

cos- = -, (13) 
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near which e l g(k) is quadratic. As it is demonstrated below, the one-magnon branches are not renormalized at H > H s 
within our approach and Eqs. ( 1 1 )— ( 12 ) reproduce the one-magnon spectrum of the model ([T]) obtained before^ij 
(remember that in our notation the distance between neighboring spins is equal to 1/2). In particular, one has from 
Eq. ( 12 ) for the field value at which the one-magnon spectrum becomes unstable 



H c = 2 J - 1 



1 

8J' 



(14) 



In contrast to b- and c- particles (or a- and j3- ones), a-particles are of the two-magnon nature (see Eqs. p]))- 
It is one of the main findings of the present paper that their spectrum derived below coincides at H > H s with 
the low-energy part of the two-magnon bound-state spectrum obtained before by using other approache d 8 * 17 * 18 ! (see 
Fig.§. 



C. Diagram technique 

We find it more convenient not to use the unitary transformation ^ in the following and to introduce four Green's 
functions 

G a (k) = -i(a k al), (15) 

G b {k) = -i{b k b\), (16) 

G c (k) = -i(c k cl), (17) 

F(k) = -i(b k cl), (18) 

F(k) = -i{c k bl), (19) 

where k = (u, k) and a k is the Fourier transform of a^t). Notice that a-particle can not transform to a single b- or 
c- particle due to the spin conservation law (a-particles carry spin 2 while b- and c- particles carry spin 1). That is 
why the Dyson equation for G a {k) has the form G a (k) = G a0 (k)(l + E a (/c)G a (/c)), where G a0 (k) — (oj — e a0 + iS)^ 1 
and e a o is defined in Eq. ^ and which gives 

G a (k) = 1 (20) 

lo - e a0 - E a (fc) + id 

In contrast, b- and c- particles can interconvert that leads to two couples of Dyson equations for G b (k), F(k) and 
G c (k), F(k). We obtain for one of them 

( G b (k) = G b0 (k) (1 + S fc (fc)G b (fc) + U(k)F(k)) , 

\ F(k) = G c0 (k) (U(k)G b (k) + S c (fc)F(fc)) , ' 1 : 

where G b o(k) = G c o(k) = (us — + iS) -1 , and Eb. c , II and II are self-energy parts. The solution of the Dyson 
equations has the form 

OM = ^^, (23) 

F(k) - W F(k) (24 ) 

V{k) = (w-E k - T, b (k) + i5)(u) - E k - E c (fe) + id) - IT(fc)Tl(fc). (25) 

It should be noted that all poles of these Green's functions lie below the real axis. This circumstance gives a useful 
rule of the diagram analysis at H > H s : if a diagram contains a contour that can be walked around while moving 
by arrows of the Green's functions, integrals over frequencies in such a diagram give zero (see, e.g., Ref!^. Using 
this rule one concludes that there are no nonzero diagrams for Sf,(fc), S c (fc) and , II(fc) at H > H s which are 
determined solely by bilinear part of the Hamiltonian (JsJ) : 

E fc (fc) = £ c (fc) = 0, (26) 
n(fc) = n* (k) = Bl = -e- tk l 2 cos % . (27) 
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FIG. 3. (Color online.) (a) Diagrams for the self-energy part E a (fc) of a-particles at H > H s . It is implied that one should 
put Green's functions of b- and c- particles Gt{p), G c (p), F(p) and F(p) defined by Eqs. (16H19I instead of double dashed 
lines (there are 10 different loop diagrams for E (£;)). Triangles represent renormalized vertices for which we obtain equations 
presented in panel (b). Bare vertices are defined by Eqs. |7f and (fsj) . 



As a result we have from Eqs. (25 1— (27 1 T>(k) = (uj — e a (k) + iS)(uj — e^(k) + id), where e Qi| g(k) are given by Eqs. ( 11 ) 
and ( p~2[ ) . Another useful rule for the diagram analysis follows from the spin conservation law. As the total spin carried 
by all incoming particles must be equal to that of all outgoing particles, one concludes immediately, for instance, that 
four-particle vertexes are zero with only one incoming or outgoing a-particle. 

It is seen from Eq. ([7| that the three-particle terms of the Hamiltonian describe decay of a-particles into two one- 
magnon particles of b or c types. They give rise to the only nonzero diagrams for £ a (fc) at H > H s all of which are 
shown in Fig. [3] 



III. ISOLATED CHAIN 



H> H s 



To calculate S a (fc) one has to derive three vertexes Fi^^q, k) for which we obtain equations shown in Fig. |3 
The solution of these equations can be tried in the form 



Tj (q, fc) = Uj(k) + Sj (k) cos g, 

T 3 (q, k) = u 3 (k) + s 3 (k) cosg + iv(k) sing, 



3 



1,2, 



(28) 
(29) 



where wi,2,3(fc), Si,2,3(fe) and v(k) are real functions. After substitution of Eqs. ( |28[ ) and (29) into equations shown in 
Fig. ^h), we obtain a set of seven linear algebraic equations for 1*1,2,3(^)1 Si^.sWand v(ky The corresponding exact 
solution is quite cumbersome for arbitrary k. However, it simplifies greatly at k = and one has at H = H s 



r 1 (q,0) = r 2 (q,0) = ~COS 



2 Q 



r 3 (q,0) = 2J-l 



1 



1 + J 



+ (1 + 2 J) cos q + i sin q. 



(30) 
(31) 



These expressions are used below for the nematic phase discussion. Substituting the general exact solution for the 
vertexes into expression for S a (fc) shown in Fig.[3]ja) we lead to the following expression for the spectrum of a-particles 
at k < 1: 



e„(k) - 2H + 2 



D« 



3J 2 (2 



-4J 

J)- 



1 + J 



16(1 + J) 2 



(32) 
(33) 



which coincides (up to a factor of 1/4 in D» due to the unit cell doubling) with the spectrum of the two-magnon bound 
states obtained before within other approache d 8 ! 17 !. Then, the condition D\\ > gives from Eq. (33 1 that the minimum 
of the t wo-m agnon bound-state spectrum is at k 



if J > 0.375 that is also in accordance with previous numerical 
fmdings^E^ and analytical results^. It is not difficult within our approach to take into account also anisotropy in 
Hamiltonian ([I]) of the form S'|S , | +1 and S , |5 , | +2 . We omit it in the present consideration for simplicity. The reader 
is referred to RefP^ for the corresponding expressions for two-magnon bound-state spectra. 
The condition e a (0) = gives from Eq. (32) the saturation field value for the isolated chain 



H q = 2 J - 1 



1 



2 + 2J 



(34) 



which is larger (if J > 1/3) than the field value (14) at which the one-magnon spectrum becomes gapless. We find 



also for the Green's function of a-particles near the pole 

G a (wwe a (k),k) 



uj — e a (k) + iS ' 
1 + 2J 
(1 + J)(2J 2 + 2J+1) 



(35) 
(36) 



These expressions are used in the subsequent calculations. 

It should be noted that the bare completely fiat spectrum of a-particles e Q o is renormalized greatly by quantum 
fluctuations (see Eq. (32)). Because this renormalization comes from processes of a-particle decay into two one- 
magnon particles, one concludes that a-particle corresponds to a state that is a superposition of the initial state with 
two neighboring flipped spins and the great number of those with two flipped spins sitting on sites which can be quite 
far from each other. It is the structure of the wave function that is used in the standard method of the bound states 
analysis 



The instability of a-particles spectrum (32 



B. H < H 3 

at H < H s signifies the transition to the nematic phase. It is not the 
aim of the present paper to discuss in detail the quadrupolar phase in the isolated chain. But in view of the quadratic 
dispersion of e a (k) at H = H s and gaps in spectra of a- and /3- particles, some results can be readily borrowed from the 
theory of ID Bose gas of real particles! 45 * 46 ! In particular, one obtains (ctjaj) = ^y/(H s — H)/D\\ using exact results 
by Lieb and Linigei^^at T = and H rj H s . Diagram analysis shows that all poles of Green's functions Gb(k) and 
G c (k) lie below the real axis at H < H s and only self-energy parts acquire some corrections (as is shown below, this 
is not the case in quasi-lD models due to the presence of the "condensate" of a-particles). Then, (bjbj) — (ctcj) = 
and one obtains from Eqs. ([3| at H w H s 

1 ' - ■ - „h 




2 -<5J) = (a}a,) = - 1 p— (37) 



that is an extension of the previous isolated chain analysis. As it is demonstrated in Fig. |4j Eq. (37) describes well 
available numerical dat d 10 * 1 ^ within the range of its validity that reads as (a^aj) -C 1. Because <C H s at J > 0.375, 
the magnetization decays quite rapidly upon the field decreasing at H « H s that is also illustrated by Fig. |4j 

Many results obtained before using other methods can be confirmed (and extended somewhat) within our approach. 
In particular, one obtains at H s» H s from an expression for the asymptotic of the correlator of densities in ID Bose- 

<SJ + „(*)S|(0)) « (SI) 2 - - ( 1 + 1 ) + Axgg^g, (38) 
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FIG. 4. Magnetization of isolated chain with J = 0.5. Circles are numerical data taken from Ref.^and the line is drawn using 
Eqs. p7b, |33l and (341. 



AnpD\\, Ai is a constant and p = (ajoj) is given by Eq. (37). Eq. (38 1 coincides in the limit 



where n — > oo, i 

H —} H a with the corresponding expression derived in RefsJ^ESEll using the bosonization technique ( J 3> 1) and the 
phenomenological theory of multipolar phases. Our discussion provides explicit expressions for the "sound velocity" u 
and p in Eq. ( 38 1 . Notice that spin-1 excitations do not contribute to Eq. ( 38 ) due to the above mentioned circumstance 



that all poles of one-magnon Green's functions lie below the real axis, whereupon all integrals over energies in the 
corresponding loop diagrams give zero at T = 0. Then, the asymptotic of the ID Bose gas "field correlator" gives for 
the nematic correlation at H rs and n — > oo 



(S+(t)S+(t)S; n (0)5 2 - +1 (0)> = (ao(t)ot(O)) 



-4, 



^/\2nTiut\ 



(39) 



where A^ is a constant, that is consistent at t = with results of R e f s LL2LL3] 



It is seen from Eqs. (38) and (39) that static longitudinal and nematic correlators decay algebraically with the 
distance at T = 0. According to the exact results of the ID Bose gas theory^ this algebraic decay changes into the 
following exponential decay at small T: (S] +n S?) e" n / r <= and <5*+S , +5 2n S7„ +1 ) « e^"/ 2r % where r c u/2nT. 
Spin-1 excitations give contributions decaying much faster at small T due to gaps in their spectra. 

To conclude our brief discussion of the isolated chain at H < H s , we confirm the long-standing expectation that 
bound states of two magnons behave in many respects like hard-core bosons. A more detailed consideration of the 
quadrupolar phase in the isolated chain is complicated by the existence of b- and c- particles. 



IV. QUASI-ONE-DIMENSIONAL MAGNETS. H > H s . 

Let us take into account the inter-chain interaction. As it is shown below and as it was found before} ^ * 18 * 21 * 22 ! quite 
a small nonfrustrating interaction between chains can destroy the nematic phase and turn H = H s into an ordinary 
quantum critical point at which the condensation takes place of one-magnon excitations. That is why we consider 
in the present paper the inter-chain interaction as a perturbation. To demonstrate the main ideas we discuss the 
simplest interaction of the form (see Fig.JlJa)) 

H' = n ^ J'llm^l^m — ^ (^l2mSi/S lm + J( im S 2 iS 2m ) , (40) 
Im Ira 

where the first sum is over the lattice sites and the second one is over the doubled unit cells shown in Fig. [TJ^b) . 
Substituting Eqs. ^ into Eq. (40) one obtains 



(41) 
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SL a {k) = 



— — »- 



V 



FIG. 5. (Color online.) Diagrams for the self-energy part of a-particles which give contributions of the second order in the 
inter-chain interaction J[ at H > H s . Triangles stand for vertexes which are calculated at J[ = and equations for which are 
presented in Fig. |3[b ). Double dashed lines represent Green's functions Gb, G c , F and F which are calculated by taking into 
account also Eq. |42|. Black dots stand for bare vertexes (431. 



where £q is a constant and 

7^4 



4 c k 



1 



E ( -J>W + 5 (4 k - J'lo) (&£&k 

J( k J b\c\a- 3 + (J( k2 + J( k j 4&_ 2 C_3 
J ik 1+ k 3 aI4 a -3a-4 + ^lk 1+ k 3 (44 & -3&-4 + 44c- 3 C-4 



2VN 



E 



v^lki 



(42) 
(43) 



ki+k 2 +k 3 =0 



1 

N 



E 

ki+k 2 +k 3 +k 4 =0 

^lkj+ks + 9^ik!+k 4 ) 44 a -3C-4 + ( 4k!+k 3 + ^J'l^+^i ) 44 a -3& 



where all vectors k have two nonzero components, y and z, (see Fig 

cos/c z ). Naturally, we will assume from now on that momenta k in Eqs. (|6j)-(|8|) are 3D vectors with the only nonzero 
component x. It is easy to show that T-L'^ and Ti!^ do not lead to renormalization of the one-magnon spectrum at 
H > H s . Then, one obtains for the spectrum of /3-particles using Eqs. (p~2]) and (42) 



Hg.ga)) and J' u 
in Eqs. I 



t(kRi„ 



(44) 

2J((cos + 



(45) 



e^(k) = H + 1 - J + Jcosk x - cos y + - (J( k - f 10 ) 

that has minima at k = (±Q, 0, 0) and k = (±Q, 7T, it) if J{ < and J[ > 0, respectively, and Q is given by Eq. ( |13[ ). 

In contrast, and H 4 renormalize the spectrum of a-particles e a (k). Considering %' as a perturbation it is easy 
to demonstrate that the only diagrams contributing to e a (k) in the second order in J[ are those presented in Fig. [5] 
Straightforward calculations show that their contribution to e a (k) has the form 



<5e a (k) = — 2Z?j_(2 + cos k y + cos k 



D 



Z/1 

2\ 



(1 + J)(1 + 3J + 3J 2 ) 
2(1 + 2J) 2 



Then, e a (k) has a minimum at k = (0,0,0) near which we obtain using Eqs. (32) and (46) 

1 



e (k) = 2ff + 2-4J- 



1 + J 



- J' 10 - 8D ± + D»k 2 x + D ± k 2 ± , 



(46) 
(47) 

(48) 



where Dy is given by Eq. (33) and k_L is the projection of k on the yz plane. It is seen from Eq. (48) that e a (k) is 
quadratic at k -C 1. Notice also that the stiffness D± of a-particles in the yz plane is of the second order in J[. 

One finds from Eqs. (45) and (48) for the field values at which spectra of j3 and a-particles become gapless (cf. 
Eqs. (Tub and pi) 



H c = 2J - 1 
H, = 2 J - 1 



1 

8J 



1 



1 



4ol 
i 



i 



j' 



2J 



(49) 
(50) 



The condition of the nematic phase existence (H s > H c ) leads from Eqs. (49) and (50) to the following inequality in 
the first order in J[: 

3J 



'^ ol < 4J(J+1)- 



(51) 



One concludes from Eq. (51) that \J[\ should be much smaller than unity in order the nematic phas e can arise. 
Eqs. (45), ( [49] ) and (50) are in accordance with previous results obtained within another approachP^^ The reader 
is referred to RefP^ for a more detailed discussion of the influence of the inter-chain interaction on properties of the 
considered system at H > H s . 
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a) 



r a (k)= 





b) 




S T Y. a (k)= _^ 



O 





FIG. 6. (Color online.) a) Equation for the four-particle vertex r a (fc) of a-particles at _ff ~ H s . A sum of an infinite number of 
diagrams which remain finite at k, q — ¥ and J( — > plays the role of the "bare" vertex "/(k, q). Some diagrams for ^{k, q) are 
also presented, b) The diagram giving the leading temperature correction to the normal self-energy part of a-particles. Same 
notation as in Fig. [3] 

V. QUASI-ONE-DIMENSIONAL MAGNETS. H < H s . 

It is convenient to discuss the quantum phase transition at H = H s in terms of Bose-condensation of a-particles 
that is similar to condensation of magnons in ordinary magnets in strong magnetic field.HH 



A. Condensation of a-particles 

According to the general scheme^^, one has to represent clq as follows at H < H s 



i<j> 



Np, 



(52) 



where <f> is an arbitrary phase and p is the "condensate" density. New terms appear in the Hamiltonian after this 
transformation. In particular, the largest terms in the ground state energy have the form S£q = — (H s — H)p+T a (Q)p 2 , 
where T a (k) is the four-particle vertex of a-particles at H = H s . An equation for T a (k) can be represented in the 
form shown in Fig. rota). Minimization of S£q gives 



H s -H 

2r a (o) ■ 



(53) 



Although the above formulas are standard,E3 the situation is more complicated here. As it is shown in Fig. [6|a), one 
has to consider an infinite set of diagrams to find T a (k) which sum j(k, q) plays the role of a "bare" vertex. These 
diagrams remain finite at k, q — > 0, J[ — > because e a ^(k) have gaps at H = H s and J[ = 0. That is why one can 
set J[ = in -f(k, q). In contrast, the integral over q diverges at J[ — >• in the second term on the right-hand side of 
the equation for T a (k) (see Fig . [6 '&)). Due to the smallness of J[, small q x are important in the integral over q. As 
a result one obtains using Eq. (351 the following equation in the leading order in J[: r a (0) = 7(0,0) — 7(0,0)Tr a (0) 
that gives 



r.(o) = 

T = 



7(0,0) 



1+ 7 (0,0)T 

z 2 r dk 



1 

T 



z 2 

0.32 - > 1, 



(54) 
(55) 



where D\\, D± and Z are given by Eqs. (33), (47) and (36), respectively, Eq. (46) has been used to take the integral 



in Eq. (55), and we omit the unity in the denominator of Eq. (54) due to the large value of T. We remind the reader 



also that the range of validity of Eqs. ( 53 )— ( 55 ) is defined by the inequality p <C 1 
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According to the general theory of the dilute Bose gas condensation liil^IHSD the- spectrum of a-particles acquires a 
linear part at small momenta k and has the form 



e„(k) = ^e a (k)(e a (k) + A), 



(56) 



D ± k 2 ± . 



where A ~ p and e Q (k) = Dnk% 

Thermal fluctuations can be taken into account quite standardly ! 44 * 48 * 4 ^ The leading contribution from them to 
S£ o comes from the Hartree-Fock diagram shown in Fig. [6jb) . As a result one obtains the following equation for the 
boundary in the H-T plane between the fully polarized and the nematic phases at T <C D± : 



H s (0)-H s (T) 



3/2 ZT a (0) C(3/2) 



(57) 



where T a (0) is given by Eq. (54), H s (0) is given by Eq. (50) and C(3/2) ~ 2.61 is the Riemann zeta-function. Notice 



that at small nonzero temperature one has to put H s (T) given by Eq. ([57]) instead of H s in Eq. ( 53 1 



B. One-magnon spectrum renormalization 



Condensation of a-particles (52) leads to terms in the Hamiltonian which renormalize the one-magnon spectrum. 



In particular, those coming from the three-particle vertexes and contributing to the bilinear part of the Hamiltonian 
|6]) have the form at H ps H s 



«4 3) - VpJ2{ e ^ T ^ °) ( c i c - q + M-q) + e i0 ri(q, 0) (c qC f _ q + &t&t q 
q 

+e-^r 3 (q, 0)6 q c_ q + e^r*(q, 0)6 qC f _ q ) , 



(58) 



where Ti(q, 0) and T3(q, 0) are given by Eqs. (30) and ( [31] ), respectively, and we neglect the inter-chain interaction in 



the sum. There are also terms of the form o q o q , c q c q , CqTq and 6 q c q coming to from four-particle vertexes a^b^ab, 
a>tfac and aJc^ac. The analysis of these vertexes similar to that carried out above for T a (k) shows that contributions 
to W2 from them are of the order of p^/D±. Straightforward calculations demonstrate that contribution from these 
terms to the spectrum renormalization is negligible. Thus, we use only Eq. (58 1 below in order to make formulas more 
compact. 



To find the one-magnon spectrum renormalization we take into account terms (58) in the Hamiltonian and introduce 
Green's functions 



P{k) 



-i{bl k b{), Q(k) = -i(d k cl), R(k) = -i(cl k bl) 



(59) 



in addition to those given by Eqs. ( 16 H( fl9|). Then, one leads to a set of four linear Dyson equations for Gb, P, R and 
F that is derived and solved in Appendix |A[ Analysis of the Green's functions denominator shows that the spectrum 
of /3-particles can be represented in the vicinity of its minimum as 



^(k) =yJel(K)-pA, 

(3J - 1)(J(2 + J(19 + 8J(J(16J(3 + 4J) - 11) - 9))) - 1) 



A=- 



128J 4 (1 + J) 2 (1 + 5J) 



(60) 
(61) 



where e^j(k) is given by Eq. (45 1, A > at J > 1/3 and A w 12/5 at J > 1. 



As it is seen from Eq. (60), the spectrum of /3-particles becomes unstable at a certain field value H c < H s that 
would signify a transition to a phase with a long-range magnetic order. While the inequality p <C 1 can hold at 
H = H Cl the present discussion is not applicable for this quantum phase transition analysis because disappearance 
of the gap in the one-magnon spectrum leads to a finite damping of a-particles at all momenta. Besides, diagrams 
contributing to 7(fc, g) contains infrared divergences at fc, q — > if the one-magnon spectrum is gapless. As a result 
the above analysis should be reconsidered at H w H c , which is out of the scope of the present paper. 

It should be noted also that one-magnon excitations acquire finite damping at H < H s stemming from loop 
diagrams. The damping, however, is small at H w H s being of the order of p (that is much smaller than the gap in 
the spectra of a- and /3- particles). 
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C. Static spin correlators and nematic order parameter 



(62) 
(63) 



One obtains from Eqs. ^ and (52) at H < H s 

(si : ) = (Sij) = 0, 
(3^) = (at) = y/pe-**, 

where _L denotes the projection on the plane perpendicular to the field direction. Then, the condensation of a- 
particles signifies the formation of the quadrupolar phase without the conventional long-range magnetic order in 
which (S^jS^j) 7^ in each (double) unit cell. This condensate should appear also "between" the neighboring unit 
cells as the doubling of the unit cell we made is only a trick. To show this we calculate the value (^2j^i(j+i)) ( see 
Fig. [ljb)), where j enumerates sites in one of the chains, for which one has from Eqs. ^ in the leading order in p 



(64) 



One obtains after simple integration using Eqs. (30), (31), (A3) and (64) in the leading order in the inter-chain 
interaction (cf. Eq. (63)) 



-ief) 



(65) 



As it is seen from Eqs. (63) and (65), the nematic order parameter (S~ S~ +1 ) = (— 1) J \/pe^ 1 ^, where j enumerates 
now spins in a chain, has an " antiferromagnetic" order: its absolute value is the same for all j whereas it s ph ase 
differs by n for two neighboring sites (j and j + 1). Such nematic ordering along chains was predicted in RefsPSU, j n 
contrast, the nematic order is "ferromagnetic" in directions transverse to chains because the minimum of e a (k) is at 
k = (0,0,0) (see Eqs. (46) and (48)). The nematic ordering does not depend on the sign of J[. Notice also that this 
finding does not depend on the way of grouping of spins into couples shown in Fig. [ljb) due to arbitrariness of <f>. 
Let us consider the static spin correlator (SjSj +n ). To find it one has to calculate (^^To+n))' (^ij^2(j+n)) ' 
and (SijS^j that can be easily done in the first order in ^fp and in the leading order in J' x using 

Eqs. 



31, pOh 



311) , ( A2 ) and ( A3 1 . The result can be represented in the form 



( S j S j+n) = VP e 



J 



l + J 



(66) 



where n > 0, that reproduces, in particular, Eqs. (63) and (65) at n — 1. In much the same way, one obtains in the 
first order in p and in the leading order in J{ 



<^„> =Pcos(f) 



J 



^ 2 



S 



j+n 



=p sin 



~2 



l + J 

J 
l + J 



5-i 



/ 2J(1 + J) 
V 1 + 2J 



(67) 



(68) 



It is seen from Eqs. (66)-(68) that all the static spin correlators decay exponentially as n — > oo. This should be 
contrasted with the case of the isolated chain in which the algebraic decay of static correlator ( 68 ) is observed! 10 * 16 ! 



As a result of the exponential decay, static spin correlators have broad peaks in quasi-lD magnets instead of Bragg 
peaks which hight rises as J increases. In particular, one obtains from Eqs. (67) and (68) up to a constant 



(5"+5_ q ) ~p 
(S^) ~p 



■ 2 

sin q 



1 



4J 2 \4J(l + J) 



cos 2 q 



1 + 2J 



cos q 



2J(1 + J) 



1 - cos q 



(69) 
(70) 



where q = (q, 0, 0). Plots of these expressions are shown in Fig. [7] It is seen that the breakdown of the translational 
symmetry of the ground state at H < H s becomes apparent in the transverse spin structure factor which has a 
periodicity in the reciprocal space equal to n/d rather than 2tt/cI. 



A small inter-chain interaction gives rise to a weak dependence of correlators ( 69 1 and ( 70 ) on components of 
momenta transverse to chains.^ 
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FIG. 7. Static spin correlator at H < H s given by Eq. ( |69| l with q = (g,0,0). The inset shows the longitudinal static spin 
correlator divided by p and given by Eq. (701. Here d is the distance between two neighboring spins in a chain. 



D. Magnetization 



Using Eqs. ^ and the solution of Eqs. (All for Gb (and the similar expression for G c ) one finds after simple 
integration for the magnetization in the leading order in p and J[ 



( ^ } ~2 1 + 2J P ' 



(71) 



Notice that in contrast to ID case discussed above, there is a nonzero contribution to the second term in Eq. ( 71 ) 
from (cjc,-) = (b]bj). 



E. Spin Green's functions 



Spin Green's functions (or generalized susceptibilities) defined as 

/>oo 

X*pM = (S^,S?) u =i / dte^([S^(t),S^0)}) 



(72) 



where a,/3 = x,y,z, can be found straightforwardly using Eqs. (|3| and (22|—(27). In particular, components of v Q(9 
transverse to the field direction are expressed via Green's functions Gb, G c , F and F. As a consequence, they have 
sharp peaks at u> = ±e Q ( g(q) corresponding to one-magnon excitations. It is seen from Eqs. (3]) that the longitudinal 
spin Green's function has a contribution containing G a {q) (apart from a smooth background originating from terms 
c^c, b% and a^a in Sf and Sf): 



Xzz(w, q) p(G a (w, q) + G a (-U), q)) w - 



2Zpe a (q) 
{u> + id) 2 - e*(q) 



(73) 



that shows sharp peaks at lj = ±e a (q) and small q. Thus, the soft mode (561 can be observed in the nematic phase 
experimentally in the longitudinal channel. 
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F. Phase transitions at H < H S (T). Symmetry consideration. 



Let us discuss the breakdown of symmetry in the transition from the fully polarized phase to the nematic one at 
H = H S (T). The symmetry of the Hamiltonian is 0(2). Symmetry operations which do not change the nematic 
order parameter (S~ S~~ +1 ) = (— l)i ^fpe^ 1 ^ include a rotation by it as well as rotations by —<f> and — + 7r accompanied 
by a reflection. These operations form a discrete group which is equivalent to Z 2 ® Z 2 . Then, the phase transition 
in the quadrupolar phase corresponds to the continuous symmetry 0(2)/(Z 2 <8> Z 2 ) = SO(2)/Z 2 breakdown (see, e.g., 
Ref.^1) that produces, in particular, the massless excitations (the soft mode). On the other hand one can expect in the 
considered quasi-lD system that the broken symmetry at small field is Z 2 €5 SO(2) (as in a non-collinear Heisenberg 
XY- magnet). Consequently, the discrete subgroup Z 2 <£> Z 2 should be broken upon the field decreasing at H < H S (T). 
This breakdown of symmetry can happen in one transition (which can be either of the first or of the second order) or 
in two subsequent Ising (second order) transitions corresponding to two Z 2 subgroup breaking!^ 3 -' The latter scenario 
can be realized in LiCuV04, where two phase transitions (apart from the low-field spin- flop one) are observed at 
H < H S (T)!HH33] 



VI. SUMMARY AND CONCLUSION 



To summarize, we suggest an approach for quantitative discussion of quantum phase transitions to the quadrupolar 
phase in frustrated spin systems in strong magnetic field H w H s . Quasi-lD and ID spin-i models described by 
Hamiltonian are discussed in detail. The approach we propose is based on the unit cell doubling along the chain 
direction presented in Fig.[T|b) and on representation (3]) of spins in each unit cell via three Bose-operators a, b and 
c ([2]). Bosons b and c describe spin-1 excitations which spectra represent at H > H s two parts of the one-magnon 
spectrum (see Fig. [2j Eqs. ^ and ([9])-(12)). Spectra (II]) and (46)-d48j) of the boson a carrying spin 2 coincide at 
H > H s with those of two-magnon bound states found before within other approaches. It is the main advantage of the 
suggested approach that there is the bosonic mode in the theory softening at H — H s . This circumstance makes the 
consideration of the transition to the quadrupolar phase substantially standard. In the purely ID case, we rederive 
spin correlators (38) and (39) obtained before either in the limiting case of J 1 or using the phenomcnological 



theory and extend previous discussions by Eq. (37) for magnetization that describes well existing numerical data at 
H w H s (see Fig. [4]). In the quasi-lD model with the simplest inter-chain interaction (40), we calculate at H < H s 
spectra of the one-magnon band (60) and the soft mode (56), nematic order parameter ( |63[ ) and (65), static spin 
correlators (66|-(70) and magnetization (71) which are expressed via the condensate density p given by Eqs. ( |53[) — 



(55). At T 7^ 0, p is also expressed by Eqs. (53)-(55) in which one should replace H s by H S (T) given by Eq. (57). 
All static two-spin correlators decay exponentially with the correlation length proportional to l/hi(l + 1/J). This 
decay results in broad peaks in static structure factors (see Fig. [7]). The periodicity in the reciprocal space of the 
transverse static structure factor is equal to ir/d rather than 2tt Jd. Transverse components of the dynamical spin 
susceptibilities ( 72 ) are expressed via one-magnon Green's functions and contain sharp peaks at w equal to magnon 
energies. The longitudinal component x^(w,q) apart from a smooth background contains a contribution (73) from 
Green's functions of a-particles which have sharp peaks at u> = ±e Q (q). 

We apply the proposed approach to the analysis of the model describing the quasi-lD material LiCuV04 in which 
a transition at the saturation field to (presumably) a quadrupolar phase has been observed experimentally. Details 
of the corresponding calculations are presented in Appendix [5] Predictions of our theory are in reasonable agreement 
with the recent magnetization measurements in LiCuVO^. Our finding that H s (0) — H S (T) oc T 3 / 2 , Eqs. (56), (60), 



(66)-(71) and (73) can be checked in further experiments on this material which should confirm also that the p 



lase 



observed just below H S (T) is really the quadrupolar one. 
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Appendix A: One-magnon Green's functions in the nematic phase 



Taking into account terms ( 58 ) in the Hamiltonian which appear at H < H s , one leads to a set of four linear Dyson 
equations for G b (k), P{k), R(k) and F(k) that reads (cf. Eqs. (21)) 

f G fc (fc) = G b0 (k) (1 + e i H^/pT 1 {k, O)P(fc) + e l ^r*(k,0)R(k) + B£F{k)) , 
P(k)=G b0 (k) (e-^2^r 1 (k,0)G b (fc) + e-' ; ^r;(k,0)F(fc) +i?^i?(fc)) , 
R(k) = G c0 (k) (e- l ^r 3 (k, 0)G b (k) + B k P(k) + er^^T^k, O)F(k)) , 
{ F(k) = G c0 (k) (B k G b (k) + e**2y/pTxtfr O)R(k) + e 4 ^r 3 (k, O)P(fc)) , 



where P(k) and R(k) are defined in Eqs. @, G b0 (/c) = G c0 (fc) = G b0 (-k) = G cQ (-k) = (u - E^ + id)" 1 , E^ = 
+ (J'iu: ~ J'io)/2, E k and B k arc defined in Eq. @. Solving Eqs. ((ATJ one finds, in particular, in the leading order 
in p 



P(k) = ^pe- 1 * 
R(k) = y/pe-t* 



2Ti(k,0) (£ k 2 + |B k | 



r|(k, o)s k (££ + w ) - B*r 3 (k, o)(E' k - U ) 



(cj - e Q (k) + id) (a; + e«(k) - id)(u; - e^(k) + i6)(u> + ep(k) - id) 

r|(k, o)Bg - (c 2 - £ff )r 3 (k, o) - 4^gkTi (k, o) 

(u - e Q (k) + i(5) (a; + e a (k) - id)(w - e / g(k) + id) (a; + e i g(k) - id) ' 



(A2) 
(A3) 



where we set p — in the denominators. Expressions for G b and F are cumbersome and we do not present them here. 
Another set of Dyson equations for G c {k), Q(k), R(k) and F(k) can be considered in much the same way. 



Appendix B: Application to LiCuV0 4 

We apply in this appendix the approach proposed in the main text to the particular quasi- ID compound LiCuV04. 
While exchange coupling constants inside a chain Ji = —18.5 K and J = 44 K have been extracted from e xperim ental 
data in LiCuVC>4 quite precisely, small inter-chain interaction has been determined much less accurately! 30 * 54 ! Then, 
we take into account to first approximation only the inter-chain coupling J' 2 = —4.3 K that is shown in Fig. |TJa) L2Q154] 
Such a model for LiCuVC>4 is considered in RefpH using another approach that involves numerical calculations at 
H > H s and self-consistent calculations at H < H s . Thus, it is reasonable to compare our analytical results with those 
of Ref.EH. Notice that the inter-chain interaction J 2 makes the system two-dimensional (in contrast to J[ considered 
in the main text). 

Let us start with the case of H > H s and represent the inter-chain interaction in the form 



(Bl) 



lm 



where the first sum is over the lattice sites and the second one is over the double unit cells shown in Fig. |TJb) . 
Substituting Eqs. Q into Eq. (Bl| one obtains Eq. (41 1, where now 



%4 



1 

1 

N 



E 

k!+k 2 +k 3 =0 



J 2k! C i4 a -3 + Jz-^blbla-3 + J 2k3 a\c- 2 c-3 + 4-k 3 a i fo -2^-; 



«.ti.t 



(B2) 
(B3) 



E 

ki+k 2 +k 3 +k 4 =0 



+JL 



2k 2 +k; 



^2k 1 +k 3 a l a 2 a -3 a -4 + ^2k 2 +k 4 a l c 2 a -3C-4 + ^2 kl +k 3 a l4 a -3 & -4 
6|4 & -3C-4 + ^2k 1+ k 4 at l4 a -3fe-4 + ^2k 2+ k 3 a l & 2 a -3C-4^ , 



(B4) 



and J' 2k = 2J 2 cos k z (l + e lfex ). We find from Eqs. Q and (B2 ) for the spectrum of /3-particles 



e | g(k) = H + 1 - J + Jcosk x - cos y |l - 2J 2 cosk z e~ ik:c 



24. 



(B5) 
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One has from Eq. (B5) in the second order in J 2 for the field value at which e^(k) becomes gapless 



H r = 2 J - 1 



1 

8J 



2J' 2 



4 



[1 - 8 J 2 ) 
16J 3 



(4) 



[l - 4 J 2 + 8J 4 ) 
16J 5 



(B6) 



The spectrum e (k) of a-particles can be calculated in the second order in J 2 as it is done for J[ in the main text. 
Considering diagrams shown in Fig. [5] one obtains that e a (k) has a minimum at k = (0,7r) near which we have 



e Q (k) = 2H + 2-4J - 



1 



- 4J 2 - (A 



2 (1 + 4J(1 + J(l + J) (3 + J))) 



D X = (J^ 



1 + J " z v "" 2(1 + J) 3 (1 + 2J) 2 

J(l + J)(2 J(3 + J(23 + J(38 + J(5 + 4J((J - 1) J(5 + J) — 9))))) - 5) - 1 



(B7) 



(1 + 2J)(1 + 4J(J(£-1) 



16J 5 (1 + J) 4 (1 + 2J) 2 

i)-0 



6472J 5 (1 + J)VV?- 1 + 2J 2 

4^(1 - 2J(1 + J)(1 + 2 J) + g) 

1 + 2J 2 ) (40 J 2 + J 3 (20 - 8f) - 2(3 + £) + J(3 i- < )( 4 :- £ ) I 

where £ = 4J 2 (5 + 4J(2 + J)) — 3, D± > at J > 1.16 and D ± sa 1/8J at J > 1. One finds from Eq. {B7) 

1 



(B8) 



H. = 2 J - 1 



2J' 2 



2 _(1+4J(1 + J(1 + J)(3 + J))) 



(B9) 



that gives i? s 



2 + 2J ' ' 4(1 + J) 3 (1 + 2J) 2 

47.08 T in LiCuV04 in accordance with the value of 47.1 T obtained in RefM We derive from 



Eqs. (B6) and (B9) for the binding energy of two magnons given by 2{H S — H c ) the value of 0.031J which is in good 
agreement with that of 0.030 J obtained in Ref.^. Then, we find using Eqs. (B6) and (B9) that the nematic phase 
can arise (i.e., H 8 > H c ) if | J 2 \ < 7.15 K that is in accordance with the inequality | J 2 \ < 6.97 K obtained in RefPS. 
It should be noted that small discrepancies in these values are attributed to the fact that the inter-chain interaction 
is taken into account exactly at H > H s in general equations which are solved in Ref! 2 ^ numerically whereas our 
analytical results are obtained in the second order in J 2 . 

It becomes very important at H < H s that the model we discuss is two-dimensional. General result for the 
condensate density in 2D Bose gas^ reads in our notation at H w H s as 



8tt/D^DI 



(H s - H) In 



H s - H 



where Z is given by Eq. (36). Eq. (B10) is valid at 



H s -H < JD»D± 



(B10) 



(Bll) 



One obtains y / D^D± w 0.37 T in LiCuV0 4 using Eqs. Q and ( |B8| ). We find from Eqs. ^ and (|BTo| for the 
static susceptibility 



X(H) = 



dH N 



J? 



z 2 



2 J 4n y /D ll D ± 



1 + ln 



H s -H 



(B12) 



The static susceptibility has been measured recently^! in LiCuV04. One of the main results was an observation of a 
cusp in x(H) upon entering into the nematic phase with the field increasing. Then, it was found that x{H) ~ M sa t/2H S 
in the nematic phase, where M sat is the saturation value of the magnetization per ion. This finding turns out to be 
in agreement wit h the prediction of Ref.^ according to which x(H ) ~ -54M sa t/H s . 

Although Eq. (B12) is valid quite close to H s , when inequality ( |B11 1 holds, one can estimate x(H) at if s — i? ~ 
■JDhDx using Eq. (B12). We obtain from this equation by discarding the logarithm x{H) ~ 0-^0M sat / H s , w here 
M sa t = 1/2 is implied. This result agrees well with the experiment in view of the fact that Eqs. (B10) and (B12) are 
valid at H s — H w ^D\\D± up to a constant of the order of unity. 

It should be noted that a small interaction between chains making the system three-dimensional would screen 
the logarithm in Eqs. (B10) and (B12) and stabilize the long-range nematic order at finite T. Our finding that 
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H s (0) - H S (T) oc T 3 / 2 , Eqs. d56l), d60l), (|66b— (TtTJ and dT3b can be checked in further experiments on LiCuV0 4 which 



should confirm also that the phase observed just below H s is really the quadrupolar one. 
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